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Abstract 

Developing stable and robust high-order finite difference schemes requires mathe- 
matical formalism and appropriate methods of analysis. In this work, nonlinear en- 
tropy stability is used to derive provably stable high-order finite difference methods 
with formal boundary closures for conservation laws. Particular emphasis is placed 
on the entropy stability of the compressible Navier-Stokes equations. A newly de- 
rived entropy stable weighted essentially non-oscillatory finite difference method is 
used to simulate problems with shocks and a conservative, entropy stable, narrow- 
stencil finite difference approach is used to approximate viscous terms. 
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1 Introduction 

The state of numerical solutions to nonlinear conservation laws is far from complete. 
While it is commonplace to use high-order numerical methods to calculate efficient 
and accurate solutions for smooth problems, solutions of problems with shocks are 
considerably more difficult to simulate. Solution methods for these problems are 
typically high-order adaptive [1,2] or hybrid [3] schemes or highly dissipative low- 
order methods. Many methods have been devised that attempt to balance accuracy, 
added dissipation, and efficiency. Most of these methods are designed using linear 
analysis of linearized equations that do not admit the formation of shocks and thus 
do not correctly account for the character of the underlying nonlinear problem. 
Additionally, stability proofs that rely on linear analysis are dependent on the reso- 
lution and do not guarantee stability for under-resolved regions. To overcome these 
limitations, we seek numerical methods that are based on nonlinear analysis. 

Any numerical method applied to problems that admit shocks should provably 
recover the weak solution of the conservation law upon convergence [4], It should 
be further proven that the weak solution recovered is the physically realizable en- 
tropy solution [4,5]. The second condition is an uncommon property in high-order 
methods. Recent advancements in this area for the compressible Euler equations 
facilitate incorporation of these properties into high-order formulations. Tadmor 
constructed entropy consistent second-order finite volume schemes that conserve 
discrete entropy [5,6], while LeFloch and Rode [7] extended these schemes to high- 
order periodic domains. These schemes have been made computationally tractable 
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for the Navier-Stokes equations through the work of Ismail and Roe [8]. A method- 
ology for constructing entropy stable schemes satisfying a cell entropy inequality and 
capable of simulating flows with shocks in periodic domains has been developed by 


Fjordholm et al. [9] Herein, an alternative approach is developed based on a finite- 
domain entropy stability proof, which yields entropy stable methods with formal 
boundary closures. 

In this work, we have constructed high-order entropy stable finite difference 
schemes for finite domains by first developing a formal set of conditions based on a 
generalized summation-by-parts property. The entropy consistent scheme for conser- 
vation laws developed by Tadmor [6] is extended to high-order with formal boundary 
closures. Based on this new entropy consistent scheme, we develop an entropy sta- 
ble correction for dissipative numerical methods such as weighted essentially non- 
oscillatory (WENO) for simulating problems with shocks. Additionally, we have 
derived a narrow-stencil, high-order viscous operator for approximating the viscous 
terms in a provably entropy stable manner. 

Using the methodology developed herein, we demonstrate the robustness and 
accuracy of the resulting entropy stable WENO operators using Burgers equation 
and the Euler equations. We also show how schemes developed using linear stability 
can fail in the presence of a shock by comparing the newly developed entropy stable 
schemes with the energy stable WENO scheme of Fisher et al. [10]. 

Results of the present work warrant investigation into the extension of the cur- 
rent entropy-stable numerical methods into generalized curvilinear coordinates. 

The organization of this document is as follows. The theory of entropy analysis 
for finite difference methods is detailed in Section 2. Conditions and corresponding 
methods for satisfying entropy stability on finite domains using arbitrarily high- 
order accurate finite difference methods are developed in Section 3. The application 
of these methods to Burgers equation and the compressible Euler and Navier-Stokes 
equations is illustrated in Section 4. Finally, the accuracy and robustness of the 
resulting high-order schemes are demonstrated in Section 5, and conclusions are 
discussed in Section 6. 


In this section, we introduce the theory of entropy stability and define the necessary 
finite difference nomenclature for conservation laws. 

2.1 Nonlinear Conservation Laws and the Entropy Condition 

The most general form of the one- dimensional inviscid conservation law on a bounded 
domain is the integral form, 


2 Methodology 



( 2 . 1 ) 
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where q denotes a scalar or vector of conserved variables, / is the nonlinear flux 
function, and the domain bounds have been assumed fixed. As noted by Lax [11], 
solutions satisfying the integral form in 2.1 are generalized or weak solutions of the 
conservation law and do not need to be smooth or even continuous. For smooth 
problems, the strong differential form of the conservation law can be written as 

Qt + f(q)x = 0, xe[x L ,x R ], te[0,oo). (2.2) 

The solution to 2.2, referred to as a strong solution , is smooth and unique but may 
not exist for all time if the physical solution becomes discontinuous. The strong 
solution also satisfies 2.1. 

Piecewise continuous solutions to the integral form of the conservation law must 
satisfy the strong form on either side of a discontinuity. Additionally, the Rankine 
Hugoniot relation holds across discontinuities [11], 

(2 ' 3) 

where T^ denotes the discontinuity location and denotes the propagation speed 
of the discontinuity. These characteristics are derived directly from the integral 
form. 

Note that weak solutions in general may not be unique [11,12], and that only 
the physically realizable entropy solution is of interest. This solution is described 
through a limiting process of a regularized conservation law that admits a strong 
solution for all time, q £ (x,t), satisfying [12] 

q £ t+f(q% = s(f (v) (q £ ,q £ x )) x , (2-4) 

where e > 0. The viscous term on the right side of 2.4 serves as an entropy dis- 
sipative regularization (defined below) [12], making all discontinuities theoretically 
resolvable. The entropy solution satisfies 

q(x, t) = lim q e (x, t). (2-5) 

e-s-0 

The entropy solution is so named because it satisfies the entropy condition, which 
for gas dynamics becomes a statement of the second law of thermodynamics. The 
general mathematical definition of entropy is a nonlinear scalar function, S(q), with 
a corresponding entropy flux, F(q ), defined by the differential relation [13] 

S q f q = F q . (2.6) 

The mathematical entropy is convex, meaning that the Hessian is positive definite, 

( T S qg ( >0, VC + 0, (2.7) 

and yields a one-to-one mapping from conservation variables, q, to entropy variables, 
w T = S q . Premultiplying the regularized conservation law in 2.4 by the entropy 
variables yields the entropy equation, 

S q q £ t + S q f(q% = Sf + F £ = eS q fW. (2.8) 
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The viscous terms in 2.8 can be rewritten as 

zS q fx V) = z™ T fx V) = £ { wT f {v) ) x ~ £w xf iv \ (2-9) 

and because an entropy dissipative regularization [12] requires that 

> 0, Vw. (2-10) 

then, an entropy dissipative regularization ensures that entropy is always dissipated 
by the viscous terms. Substituting the definition in 2.10 into 2.9 yields 

This relation is substituted into 2.8 to find the local entropy inequality [11], 

S$ + F*<e(w T fW) x . ( 2 . 12 ) 

The entropy condition for the conservation law is found by integrating 2.12 over 
space and taking the limit e — > 0, 

Xr 

jf/sd* + f|«<0 

XL 

If the weak solution satisfies 2.13, then it is the entropy solution consistent with the 
definition in 2.5. It is important to note that the mathematical entropy has the op- 
posite sign from thermodynamic entropy in gas dynamics. Thus, the mathematical 
entropy across a shock decreases instead of increases. This nomenclature is used 
consistently throughout this document. 

2.2 Entropy Analysis 

The application of continuous entropy analysis was used above in the derivation 
of the entropy condition. Some additional formal definitions are useful to further 
specify the mathematical characteristics of the entropy. 

As stated in Section 2.1, the mathematical entropy is a nonlinear function of the 
conservation variables, S(q), with a corresponding nonlinear entropy flux, F(q). A 
set of entropy variables, w, with a one-to-one mapping to the conservation variables, 
q, is defined based on this entropy. The entropy variables have some remarkable 
properties. Because of the one-to-one mapping, the conservative variables can be 
written as a function of the entropy variables, q(w). Thus, the strong form of the 
conservation law can be rewritten as, 

q W Wt + fqQwWx = 0. (2-14) 

Both q w and f w = f q q w are symmetric matrices [13], so the conservation law is a 
symmetric hyperbolic system when written in terms of entropy variables. 


(2.13) 


( 2 . 11 ) 
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In the development of the entropy condition, we assumed that the regularization 
terms were entropy dissipative, satisfying 

> 0, Vw. 

Following Hauke et al. [14] this is shown by casting the viscous flux in a quasi-linear 
form, 

f (v) = c{q,x)q x , (2.15) 

where c(q, x ) may be a scalar constant, a variable matrix that depends on spatial 
location, or nonlinear in the conservation variable. An entropy function can be 
chosen such that the viscous coefficients are symmetric and positive semi-definite, 

f {v) =cw x , c = c(q,x)q w , ( T c( T > 0,VC (2.16) 


It is clear that if this transformation holds, then 

w xf ( ' V ' > = w x cvj x > 0, \/w, (2-17) 


and the viscous regularization is indeed entropy dissipative. Other authors [6,13,14] 
have noted that an entropy that makes the hyperbolic matrices symmetric will not 
necessarily make the diffusive coefficient matrix symmetric. Therefore, the space 
of possible entropy functions for a parabolic or incompletely parabolic system is 
reduced compared to the hyperbolic problem. This consideration is important when 
defining the entropy condition for a given problem. Herein we restrict our definition 
of entropy stability to entropy functions that satisfy a specific viscous regularization 
even when evaluating the stability for problems in the limit of zero viscosity. 

For nonzero viscosity, the continuous entropy decay rate is found by substituting 
2.16 into 2.8 and integrating over space, 


dt 


XR 


j S dx = 

ew T f W - F 


XR 
J X L 



XL 


XL 


(2.18) 


The last integral term is positive semi-definite and thus the entropy will only increase 
in the domain through the boundaries. The goal of the numerical methods designed 
in this paper will be to mimic 2.18 at the semi-discrete level. 

Harten [13] describes that the symmetry of the matrices, q w and f w , indicates 
that the conservation variables, q, and flux, /, are Jacobians of scalar functions with 
respect to the entropy variables, 

q T = <Pw, f T = i’w, (2.19) 


where the nonlinear function, <p, is called the potential and ■0 is called the potential 
flux [6]. These nonlinear functions satisfy 

ip = w T q — S, ip = w T f — F. (2.20) 


Just as the entropy function is convex with respect to the conservative variables 
( S qq is positive definite), the potential function is convex with respect to the entropy 
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variables. We note that the one-to-one mapping admits an alternate form of the 
flux based on the entropy variables, g(w) = f(q ). 

The entropy and corresponding entropy flux are often referred to as an entropy- 
entropy flux pair, ( S,F ). Similarly, the potential and the corresponding potential 
flux are referred to as a potential- potential flux pair, (ip, ip) [6]. The symmetry 
properties and the definition of the potential flux are used in the stability analyses 
in the rest of this work. 

Entropy analysis is valid for nonlinear equations and discontinuous solutions. 
It is therefore more generally applicable than linear energy analysis and gives a 
stronger stability estimate. We now turn our attention to the mathematical for- 
malism required in spatial discretizations in order to mimic the continuous entropy 
properties at the semi-discrete level. 

2.3 Spatial Discretization 

Most finite difference approximations rely on a uniform discretization of the domain 
in each direction. Typically this uniform discretization is conducted in a compu- 
tational space, and then a transformation to a nonuniform physical space permits 
greater flexibility in approximating solutions with varying scales. In this work, 
we limit our attention to Cartesian domains and extend the results to curvilinear 
multi-block domains in the future. 

An important element in the approach taken here is the use of complementary 
grids. These grids allow the finite difference operations to be written as simple flux 
differences, analogous to the approach of the finite volume method. In a previous 
paper [15], we showed that this telescopic flux difference form yields a generalized 
summation-by-parts (SBP) property that is used to show that the weak solution to 

2.1 is recovered when the solution converges. 

2.3.1 Complementary Grids 

The domain hi = [xl, x r ] is divided into (IV— 1) uniform segments with N equispaced 
endpoints denoted by x, 

rj~i 2 1 

X = (xi,x 2 , ■ • • ,X N ) , Xi = x L + _ - (x R - x L ) , 1 = 1,2, ...N. (2.21) 

Since the approximate solution is constructed at these points, they are referred to as 
solution points. It is useful to define a set of intermediate points prescribing bound- 
ing control volumes about each solution point. These ( N + 1) points are referred 
to as flux points as they are similar in nature to the control volume edges employed 
in the finite volume method. The distribution of the flux points depends on the 
discretization operator. In standard second-order finite difference methods, the flux 
points are located half way between adjacent solution points. For higher-order fi- 
nite differences, the flux points are located half way between solution points in the 
domain interior, but the spacing between flux points abruptly becomes nonuniform 
as the boundaries are approached to satisfy the summation-by-parts condition (ex- 
plained below). The spacing between the flux points is incorporated into the finite 
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difference operator using the norm, V . where the diagonal elements of V are equal 
to the spacing between flux points, 


x = (xo, x ii ■ ■ ■ x n) T , xo = xi, xn = xn, 
Xi Xi — i ^ 1? 2, ... , N. 

In operator notation, this is equivalent to 


Ax = VI 


where 

1 = (1,1,...,1) T , 


is a vector with N elements and 


/ -1 1 0 

0 -1 1 


0 0 0 

\ o o o 


0 0 0 \ 

0 0 0 

■■•00. 

-1 1 0 

0 — 11 / 


( 2 . 22 ) 


(2.23) 


(2.24) 


is an N x ( N + 1) matrix that calculates the undivided difference of the two adjacent 
flux points evaluated at the solution point. Note that in 2.22, the first and last flux 
points are the same as the first and last solution points. This consistency is needed 
to define unique operators. The discretization is illustrated in Figure 1. 


fi-i 

Ui-l 


Xo 


+ » + 


Xi X 2 Xs x 4 

X\ X2 X$ X4 


fi fi+1 fi -\- 2 

Ui Ui+1 u i+ 2 

f ^ X l _x t _x t X I xU 

Xi XN-i xn- 3 xn- 2 xn- 1 1 

Xi Xi + 1 XN -3 Xjsr-2 Xn- 1 Xn 


Figure 1. The one-dimensional discretization for finite differences is illustrated. 
Solution points are denoted by • and flux points are denoted by x . 


The approximate solution on the grid is denoted by 

u(f) = (ui{t),u 2 {t ), . . . , u N (t)) T , Ui(t) = u h (xi,t) i = 1,2, . . . ,N, (2.25) 

where q(x,t ) has been reserved to denote exact solutions, and Uh(x,t) denotes ap- 
proximate solutions. Quantities located at flux points are denoted with an overbar. 

Finite difference methods approximately satisfy the governing equation at the 
solution points, x. Difference methods based solely on the solutions points are 
considered next. These methods are then manipulated into simple flux differenc- 
ing methods that use interpolated data at intermediate flux points. Recasting the 
methods in this way facilitates implementation as well as conservation properties. 


2.3.2 First Derivative Approximation 

We utilize first derivative approximations that satisfy the summation-by-parts (SBP) 
condition, which means that the derivative approximation mimics integration-by- 
parts, 

XR XR 

J <jm x dx = (f>u\l* - J 4> x udx. (2.26) 

XL XL 

This mimetic property is achieved by constructing the first derivative approximation, 
T>(f), with an operator in the form 


V = P~ 1 Q, V = V T , C T VC>0, C/0, 
Q t = B-Q, B = diag (—1, 0, . . . , 0, 1) . 


The function of the diagonal matrix V is to incorporate the local grid spacing into 
the derivative definition. The nearly skew-symmetric matrix, Q, is an undivided 
differencing operator where all rows sum to zero and the first and last column 
sum to —1 and 1, respectively. The difference operator, V, approximates the first 
derivative as 

0z(x) = V(j) + T P 2pp, (2.28) 

where Tp2pp is the truncation error of the approximation. The nomenclature 2 p 
refers to the interior accuracy and p refers to the accuracy at the left and right 
boundaries. Note that this truncation error ( T P 2 PP ) is optimal if V is constrained to 
the diagonal norm case on a uniform grid and still satisfies the SBP property [16]. 
Integration in the approximation space is conducted using an inner product with 
the integration weights contained in the norm V, 

XR 

J (j)u x dx « <f) T VT>u, 4> = (0(xi),0(x 2 ), ■ ■ ■ A(xn)) T ■ (2.29) 

XL 

Using the definition in 2.27, the SBP property is demonstrated, 

(p T VV 1 Qu = 4> T (B - Q t ) u = 4>nun — 4>iui - 4 > t V t V u. (2.30) 

The specific operator used in this work is shown in Appendix A.l. 


2.3.3 Variable Coefficient Second Derivative Approximation 

The viscous approximations for regularized conservation laws, written in general as 

(d(x)v x (x.)) x = £> 2 (tf)v + T^ 2 l p , (2.31) 

must also satisfy the SBP condition. Integration by parts yields 

XR Xr 

J 4>(^Vx) x dx = <j>&v x \lf - J 4>x®v x dx. (2.32) 

X L X L 
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It is trivial to show that two applications of the first derivative operator satisfy 
the SBP condition. In practice, this is not advisable, as the approximation using 
two first derivative operations requires a much wider stencil (is less efficient), is less 
accurate, and leads to only neutrally stable approximations [17-19]. Instead the 
viscous operator is defined as 

V 2 {0)=V- l {-M{#)+B[0\V), = [$] = diag (#(x)) , 

C T M(tf)C>0, C T [^]C > o, VC- 

The inner product with the 'P-norrn yields 

4 b T VV + B[d]V) v = 4> t B[$]V u - 4) T M('d)v. (2.34) 

It is clear that the boundary terms are mimicked from the continuous case, but 
based on the definition in 2.33 it is unclear that 

XR 

j 4>x^v x dx « </> r .M(i?)v. 

XL 

The easiest way to show this is to restrict the definition in 2.33 to what Mattsson [19] 
calls compatible operators, 

M{ti) = V T V[d}V + ft(tf), K(d) = 'R{V) T , C T ^)C > 0, VC, (2.35) 

where 7£(i?) is called the remainder matrix and scales with the grid spacing at 
the same order as the truncation error of the second derivative operator. This 
shows the relationship between applying two first derivative operators (called a 
wide stencil approximation) and directly applying the second derivative. lZ{d) is 
chosen such that the wider stencil terms in V 1 VT> vanish, but the approximation 
of the second derivative still holds. This is detailed below. The SBP condition is 
then demonstrated by 

4) T VV 2 (V)v = 4> T B[d]Vv - 4> t V t V[V]Vv - 4) T 1Z(^)w. (2.36) 

The last term is small and decreases quickly with increasing resolution. Therefore, 
the SBP property for the viscous term mimics integration by parts at the designed 
order of accuracy. 

To calculate the remainder matrix, we alter slightly the approach taken by Matts- 
son [19], 

n r 

n$) = Y, C T [^]C>0, VC, (2.37) 

k = 1 

where A/& are rectangular operators for higher derivatives and [$]*, are diagonal 
matrices with a convex combination of the variable coefficients. The major difference 
between this approach and the approach originally proposed by Mattsson is the use 
of rectangular matrices, which simplifies the derivation somewhat. The derivation 
of the matrices for fourth-order is shown in Appendix A. 2. 
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2.3.4 Telescopic Flux Form 

All flux gradient operations in this work are recast into a telescopic flux form, 

f x (u) = V~ 1 Ai + T P 2 PP . (2.38) 

This demonstrates the utility of the complementary grid as described in Section 
2.3.1. We show in a previous paper [15] that various conservative and accurate flux 
gradients can be constructed in this manner with formal boundary closures based 
entirely on the SBP operator, Q. This is extended to a general formula for entropy 
consistent finite differences in Section 3.4. When Q satisfies the SBP condition, the 
endpoint fluxes are always consistent, 

fo = f(u i), In = f{u N ). (2.39) 


This telescopic flux form admits a generalized SBP property. The typical SBP 
operator in 2.27 transfers the action of the discrete derivative onto a test function 
with an equivalent order of approximation. The telescopic flux form combined with 
the flux consistency condition results in a more generalized relation, 


tfW-'M = - A)f = f(u N )</> N - f(u i)<)>i - 0 Af, 


(2.40) 


where 



( 0 

-1 

0 

0 

0 

0 

\ 


( - 1 

0 

0 

0 

0 

0 

\ 


0 

1 

-1 

0 

0 

0 



0 

0 

0 

0 

0 

0 


A = 

0 

0 



0 

0 


, B = 

0 

0 



0 

0 



0 

0 

0 

1 

-1 

0 



0 

0 

0 

0 

0 

0 



V o 

0 

0 

0 

1 

0 

J 


V o 

0 

0 

0 

0 

1 

/ 


and 


^0 T A = 0j + o(iv- 1 ) 


This is equivalent to the commonly used explanation of summation-by-parts in in- 
dicial form, 


N N - 1 

^2 <t>i ( fi - fi- 1) = f(u N )(t>N - f{ui)(/>l -J2fi (<A+! - <t>i) ■ 
i= 1 i= 1 


(2.41) 


The action of the derivative is still moved onto the test function but at first order 
accuracy. This generalized property is important for satisfying the weak form, as 
shown in Section 2.4. 

The variable coefficient viscous operators described in Section 2.3.3 can be con- 
structed to satisfy 

(0v x (*))* ~ V- 1 + B[0\D) v = V~ l M [v) . (2.42) 


This is shown by construction in Appendix A. 2.1. 
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Remark. The flux form is not commonly used in finite difference implementa- 
tions, except by those that use WENO or hybrid WENO with central differencing. 
It is unnecessary to implement a scheme in this way, but it is important to be able 
to show that a scheme can be cast in this form. This was a critical step toward 
proving that different finite difference forms can be constructed to satisfy the Lax 
Wendroff theorem [15]. Furthermore, it is shown in Section 4.2.3 that sometimes the 
flux form is necessary to describe a difference operator when no simple differential 
form exists. 


2.3.5 Semi-Discretization 

The finite difference approximations described above are used to change the system 
of partial differential equations into a set of coupled ordinary differential equations. 
The result is the semi-discrete equation, 

ut + V^Ai = V~ 1 At v) +V- 1 g b , (2.43) 

where gj, contains the enforcement of boundary conditions, detailed in our previous 
work [15] and other references [20-24]. The algorithmic form, 

(ft 1 - /£>,) - (/, - ;-i) + 

is used in all applications in this work and illustrates how the complementary grids 
interact. Gradients are specified only using simple differences of interpolated fluxes 
at the flux points. 

2.4 Satisfying the Weak Form 

We showed in a previous paper [15] how the telescopic flux form described above 
can be used to guarantee that the weak solution is recovered when the solution 
converges. Two additional conditions are required [25]: All fluxes in f must have 
compact support, 

fj = fj {uj-i+i,Uj-e+2, ■■■, Uj+e-i,Uj+i) , j =0,1,..., iV, (2.45) 

where i > 0 is a finite integer independent of N ; and f must be consistent with the 
flux in the Rankine Hugoniot relation, 


( u i)t 


( 9b) i 

, i = l,2,...,N, (2.44) 


fj (q,q, ■ ■ ■ ,q,q) = /(<?)■ (2.46) 

In other words, as the resolution increases, the flux stencil width must remain con- 
stant, and the functional form of the flux must match the functional form of the 
conservation law flux. 
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2.5 Temporal Integration 

In all simulations used herein, the low-storage, five-stage, fourth-order, Runge Kutta 
scheme of Carpenter and Kennedy [26] is used to advance the semi-discretization in 
time. It is noted that this scheme does not satisfy the strong stability preserving 
(SSP) property [27]. In the problems of interest for the current work, schemes that 
did satisfy the SSP property exhibited the same robustness and stability as the low- 
storage scheme, so the extra cost associated with these methods was unnecessary. 
However, SSP Runge Kutta time integration and an appropriate limit on the time 
step are required for the formal proof to hold in time. 

3 Entropy Stable Finite Differences 

In this work, we seek spatial flux divergence approximations using the finite dif- 
ference method that mimic the continuous entropy condition at the semi-discrete 
level. This is accomplished through a semi-discrete entropy analysis to determine 
the conditions on the telescopic flux form required to satisfy the mimetic property. 

3.1 Semi-Discrete Entropy Analysis 

The goal of the semi-discrete entropy analysis is to show that the semi-discrete 
entropy mimics a continuous entropy estimate such as that given in equation 2.18. 
The nonlinear analysis begins by contracting the entropy variables w 2 with the semi- 
discrete equation 2.43. The resulting global equation that governs the semi-discrete 
decay of entropy is given by 

+ w r Af = w r Af (t ^ + w T g^, (3-1) 

where 

w = (w(ui) T ,w(u 2 ) T , ■ ■ • , w(u n ) t ) , 

is the vector of entropy variables. The time derivative is easily manipulated into a 
mimetic form by using a diagonal norm SBP operator (which commutes with any 
arbitrary diagonal matrix) and the pointwise definition of entropy 

w? (ui) t = (Si) t , Mi. 


The resulting expression is 

w t Ru t = 1 T, PS f 

Using this result, equation 3.1 is recast as 

(3.2) 

This equation and mimetic arguments for semi-discrete entropy consistency will 

— — ('ll) 

yield sufficient conditions for f, f , and g 5 . The derivation of entropy consistent 
inviscid and viscous terms, however is considerably more involved than that required 
for the time term. 
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3.2 Inviscid Flux Conditions 
3.2.1 Entropy Consistent Fluxes 

By definition, an entropy-consistent inviscid semi-discretization satisfies the expres- 
sion 

w T Vu t + F(u n ) - F(m) = + F(un) ~ F(m) = w 1 g b . (3.3) 

Comparing this expression with inviscid portion of equation 3.2 immediately reveals 
that the inviscid flux terms are entropy consistent if they satisfy 

w r Af = F(u n ) — F(ui) = 1 t AF. (3.4) 

It is difficult to enforce the global entropy consistency condition in 3.4 at the flux 
points. To circumvent this problem, Tadmor [6] developed a more restrictive condi- 
tion that yields a local condition on the flux for the global entropy consistency. We 
generalize this condition in the current work for higher approximation orders and 
finite domains. 

Substituting the definition for generalized summation-by-parts in Section 2.3.4, 
A = B — A, into the global entropy consistency condition in 3.4 yields 

w T 6f - w T Af - i t bf + i t af = w T m - 1 t bf - w T Af = o. ( 3 . 5 ) 

The boundary terms in 3.5 can be reorganized as 

w 1 Bf - 1 t BF = (w%f N - F n ) - (wjfi - Fi) = ip N - -01 = ^ Bl, 

where ifti and i/jn represent the potential flux defined in 2.20, and ip is a vector of 
potential fluxes, discussed below. Defining [f] as a diagonal (N + 1) x ( N + 1) matrix 
containing the elements of f, 3.4 further simplifies to 

( tp T B - w T A[f]) 1 = 0. 

~T ~T ~ - 

Substituting the equality tp 131 = ip A1 into the left side of the equation yields 

(ip T A — w T A[f]^ 1 = 0. (3.6) 

This is satisfied by the vector sufficient condition, 

Xp T A = w r A[f], ^N = ipN i (3.7) 

and subsequently by the local sufficient conditions 

(u’i+i ~ Wi) T fi = A+i ~ V>i, * = 1 , 2 ,..., IV- 1. (3.8) 

A flux that satisfies this condition given in equation 3.8 is denoted f ‘ . 

Note that Tadmor arrives at the condition ip = ip because of an assumption 
on the form of F [6]. In the generalized condition derived in equation 3.7, it is 
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unnecessary to define ip in the domain interior. Indeed, it is not even unique because 
of the arbitrary assumptions used to relate the matrices in the domain’s interior. 
We shall see next that this generality is important for high-order methods, because 
the consistent entropy flux does not satisfy Tadmor’s original form [7]. 

The following theorems are two important contribution of this work. They prove 
that high-order entropy consistent fluxes can be constructed from linear combina- 
tions of two-point entropy consistent fluxes. This results follows immediately from 
the structural properties of diagonal norm SBP operators, because they all admit 
the generalized SBP given in section 2.3.4. 

Theorem 3.1. A two-point entropy consistent flux can be extended to high order 
with formal boundary closures using the form 

N i 

fi S) = ‘ 2( htk).fs (ue, U k ) , 1 < i < N - 1, (3.9) 

k=i- 1-1 1=1 

when the two-point non- dissipative function from Tadmor [6] is used 
i 

Is ( u k , ufl) = J g ( w(u k ) + f ( w(u t ) - w(u k ))) df , g(w(u )) = f(u). (3.10) 

o 

The coefficient q( k ,£) corresponds to the ( k ,£ ) row and column in Q, respectively. 
Proof. To show the accuracy of approximation, we express the flux difference as 
N i N i— 1 

ip-fn= e e 2q(e,k)fs (ue,u k ) - EE 2 q(e, k )fs {u£,u k ) , 2 < i < N — 1. 

k=i -\- 1 1=1 k—i 1=1 


The stencils of each flux contain considerable overlap, which is apparent if we rewrite 
the difference as 


N i-1 


N 


fz S) ~fi S l= ^ 2 q(i,k)fs(ui,u k ) 

k=i-\-l 1=1 k=i+ 1 

N i—1 i—1 

- EE 2 <?(£,a -)fs (ue, u k ) - ^2 2 q(£,i)fs (ue, uf ) , 


£=\ 


k=i + 1 1=1 

N i-1 

: X! 2( l(hk)fs (Ui,u k ) - ^2 2 q(£,i)fs (ue, uf) 
k=i -\- 1 t=l 

2 < i < N - 1. 


(3.11) 


Noting that q^j) = b ^ - q^ and that fs(u k ,ue) = fs(u£,u k ), 3.11 is rewritten 
as 

N 

j(S) _jKS) = -£ 2 q(i>j) f s ( Ui , Uj ), 2 < * < N - 1. (3.12) 

3 = 1 


15 


By the same argument, the left boundary flux difference is 


N N 

f[ 5) - fo S) = ^2 2 Q(i,k)fs (ui,u k ) - /Oi) = 22 2q ( 1 ,k)fs (U!,u k ) , 

k = 2 fc= 1 

and the right boundary difference is 


N- 1 N 

In'* “ In-i = f( u n) ~ 22 2 q {t,N)fs (ui, u N ) = 22 2 Q(N,k)fs (u N , u k ) . 

e=i fc=i 

Combining the expressions for the boundary points with 3.12, the flux difference at 
all solution points is expressed as 

N 

fi S) - fi- 1 = *22 2( l(i,j)fs( u iiUj), 1 <i<N. (3.13) 

3 = 1 

This form facilitates an analysis by Taylor series at every solution point. The accu- 
racy constraints of Q satisfy 


E 


Qik(k - if 

l\ 


6x 


he, 


e = o,i,...,d, 


(3.14) 


where d = 2p in the domain interior and d = p at the boundaries. To proceed 
further, we need to rewrite the flux function in 3.10 as a Taylor series expansion 
about £ = 0 [5], 


9 ( ) + £ (w(uj) - w(m ))) = 222 


1 d k g 


'4—' k\ dw k 
k = o 


k /- k 


(Wj - Wi ) '£ 


Wi 


The integration now proceeds simply, 

7 / \ f ST' 1 d k 9 

k= 0 


0 


E 

k = 0 

oo 

E 

k = 0 


1 d k g 


k\ dw k 

1 d k g 
k\ dw k 


( Wj - Wi) k i k d£ 
l 

(wj - Wi ) k J £ fc d£ 

o 

1 


(3.15) 


(Wj - WiY 


k + 1’ 


and simplifies to 


fs(ui,Uj ) = g(wi) + 22 


d k g 


— ( [k + 1)! dw k 


(Wj - Wif 


Wi 


(3.16) 
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We now perform a second Taylor expansion about Wi, 


f S {ui,Uj) = g( Wi ) + ^2 


1 d k g 


E 

w i ll= 1 


(< 5 x)^(j — d f 'w 


J (A: + 1 )! ch/; fc 

This expression for the flux is substituted into 3.13, 

N 


t\ 


dx £ 


-i k 


Xi 


(3.17) 


fT-i H = E 2 «> 

i=i 


+ y 


1 


— j (fc + 1 )! dw k 


E 

U=1 


( Sxf(j — if d l w 


e\ 


dx £ 


-i k 


Xi 


The first term vanishes because all rows of Q sum to zero. We then split k = 1 from 
the remaining terms in the sum, 


N OO rs 

AS) _ AS) d 9 

J i J i - 1 /_ ^ Z_> g w 

j = i i=i 

N oo 


( Sx) £ qij(j — iy d e w 


a 


dx l 


_l_ ^ <9 fc g 


(A: + 1 )! < 9 u> fc 

1=1 fc =2 v ’ 

Using 3.14, the first term is rewritten as 



'y> (fe)^(j - i) £ 


Wi 

u=i 

rr*_ 


(3.18) 


N oo 

EE 

j=i e=i 


dg 

dw 


i\ 


i) £ < 9 £ ie 

dg 

dw 

dx e 

T . 

u-'Z 

f) r 

Wi UJb 


V (i)(i) + ° 


d+1 


which will yield a design-order accurate approximation of g x = f x at x,;. We now 
must show that the remaining sum term in 3.18 is design-order small. To handle 
the series raised to a power, we recursively apply the identity, 


y akxk I ( y b k 


x k \ = 


\k = 1 


\k = 1 


y ctfx 

Ui=l 


7i 


y ^ 


OO OO 


X ^ = 


y y 


,(<i+«0 


U 2 =l 


ti=ie 2 =i 


such that the expansion to the power k can be rewritten as 


E 

U= l 


Using 


( Sxf(j — if d^w 


t\ 


dx l 


-i k 


OO oo 


yy | d (rn w\ (6x)^^ m (j — i)^^ n 
i \= 1 = l \m=l 


cU/” 


(E- 


y > 1 , VA: > 2 , 


m= 1 


we can apply the accuracy constraints 3.14 to the second sum term in 3.18, 


N OO 


yy yy 2qij 


l=i fc=2 


(A: + 1 )! dw k 


y~-v ( 5x) e (j — i) e d e w 
v i le=i 


i\ 


dx e 


= O [(5x) d+1 ) . 
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Thus we have shown that using 3.9, 


AS) _ AS) 
Jj n - 1 

v m) 


f x + 0 (( 5x ) d ) . 


(3.19) 

□ 


Theorem 3.2. A two-point high- order entropy consistent flux satisfying 3.8 with 
formal boundary closures can be constructed using 3.9, 

N i 

fi S) = ( U A Uk), 1 < i < N - 1, 

k=i -\- 1 £= 1 

where fs(ue,Uk) is any two-point non-dissipative function that satisfies the entropy 
consistency condition 


(we - w k ) T f s (u e , u k ) = ife~ 'f’k- (3.20) 

The high-order entropy consistent flux satisfies an additional local entropy consis- 
tency property, 

WV^Af^ =r~ 1 AF = F x (u) + T d , (3.21) 

or equivalently, 


wf(ff ) -fl S . ) 1 ) = (Fi-F i _ i), l<i<N, (3.22) 


where 


N i 

Fi = [(wt + Wk) 1 fs(ue,u k ) - (i>e + 4>k)\ , 1 <i < N - 1. (3.23) 

k=i -\- 1 £=1 


Proof. Using 3.13, the inner product of the entropy variables with the flux difference 
can be expressed as 


N N 


W T Af ‘ 2q (i,j) W f fs ( u i > u 3 ) 

i= 1 j = 1 
N N 

= ££ (q(i,j) + \i,j) - 9(j.<)) wffs(Ui,Uj). 
i= i j = i 


(3.24) 


Using the structure of B and recognizing that the summation indices are arbitrary, 
and thus, q(jj) w i = Q{i,j) w ji we rewrite 3.24 as 


N N 

W 7 Af (5) = wJjf N -wjfi + EE Q(i,j)( w i - Wj) T fs(ui,Uj). (3.25) 

i = 1 3 = 1 
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Tadmor [5] proved that 


(m - w k ) T / g ( w{u k ) + i ( w(u e ) - w(u k ))) d £ = ip e - ip k , 


(3.26) 


but more generally we can show that any flux satisfying 3.20 can be used to simplify 
3.25 to 


w 


N N 

’ r Af ( ‘ S) = w NfN ~ wjfi + EE««>c*-*> 

i=l j = 1 

N N N N 

WnIn -wjfi + 

2=1 j= 1 2=1 j= 1 

N N 

Fn + ipN - Fi - ipi - q(i,j)ipj = F n - Fi , 

2=1 J = 1 


where we have used w T f = F + tp. To the authors’ knowledge, this proof of entropy 
consistency for high-order finite difference methods on bounded domains is new. 

The accuracy proof of the local entropy consistency property 3.21 follows imme- 
diately from the accuracy proof of the flux difference 3.19, 




V, 


= w. 


(*)(*) 


f (fx(ui) + O ((fe) d )) = F x (ui) + O ( (fey 

1 <i< N. 


(3.27) 


To show that 3.9 satisfies the local entropy consistency property 3.22, we start by 
writing the entropy flux difference, 


N i 

Fi - Fi_ 1 = "^2 q (W) + w k) T J s (w, u k ) - (tpe + i/} k )] 

k=i-\-l £=1 
N 2-1 

-EE Q{i,k) [( w £ + u ’ k ) T fs (u e , u k ) - + ip k )\ , 

k—i i—\ 


2 < i < N - 1. 


(3.28) 


Following the procedure in 3.11 to get 3.13, we similarly get the simplified form of 
the entropy flux difference 

N 

Fi — ^ [(u>i + Wj) T f s (ui,Uj) — (ipi + ipj)] , 1 < i < N. (3.29) 

3 = 1 

To show the consistency, we premultiply 3.13 by wf , 

N 

W I ( fi S) - fi-l) = 2 <l(i,j) w i fs{ui, Uj ), 1 < i < N. (3.30) 

3 = 1 
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Adding and subtracting 


N 

^2q(i,j)wJ fs(ui,uj) 

3 = 1 

and using 3.20, we can rewrite the difference as 

N 

wf ( f- s) - = Y ?(ij) [(“>» + Wj) T fs(ui,Uj) + (Wi - Wj) T fs(u uUj )] 

3 = 1 
N 

= X! [(™» + w j) T fs( u i, Uj ) + (V’i - Y,)] • 

i=i 

We now simply add 

N 

i=i 

to get the flux difference formula consistent with 3.29, 

N 

wj (]\ S) - jjfj) = X] [( w < + w o)fs (ui,uj) - (V'i + Yb')] = Fi- Fi- 1 , 

j'=l 

1 < i < N. 


w 

Using Theorems 3.1 and 3.2, we are guaranteed that the extension of the two- 
point flux, 3.9, is a high-order accurate entropy consistent discretization of the 
conservation law. 

Remark. The entropy consistency proof is satisfied for all two-point fluxes that 
satisfy 3.20. The accuracy proof has only been proven for fluxes in the integral 
form 3.10. We are at this point unable to show that any flux satisfying 3.20 will be 
design-order accurate, so such fluxes should be validated for accuracy independent 
of Theorem 3.1. 

3.2.2 Entropy Stability 

For entropy stability an analogous condition to 3.3 is 

w T Vu t + F(u n ) - F(ui) < w T g b . (3.31) 

We again substitute the semi-discrete conservation law into 3.31 and show that the 
entropy stable inviscid fluxes must satisfy 

w T Af > 1 t AF. (3.32) 

Using the result for the entropy consistent flux in 3.4, this condition can be rewritten 
as 

w 1 Af > w T Af f ' S \ 
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Substituting the generalized summation-by-parts property, 


w t (B - A)(f (5) - f) = w T A(f - f (S) ) = w T A([f] - [f {5) ])l < 0, (3.33) 

the sufficient local condition for entropy stability is 

w T A[f] < w T A[f (S) ], (3.34) 

or in indicial form, 


{w i+l - Wif (fi - f\ S) ) <0, * = 1,2,..., iV — 1 . (3.35) 

The development of this condition relied heavily on the formalism introduced in 
Section 2.3. The generalized SBP property made it possible to extend the entropy 
stability condition to finite domains. 

Remark. The entropy stability condition is based on the global property 3.32. It 
has not been proved to provide a pointwise entropy stability property at this point. 

The entropy stability condition in 3.35 can be used to find entropy stable con- 
ditions for any type of telescopic flux operator. It only informs the upper limit of 
dissipation from a flux operator such that the entropy condition is satisfied. It does 
not inform how to add sufficient dissipation such that a non-oscillatory solution is 
obtained. In this work, we use the WENO finite difference method to construct 
dissipation and use the condition 3.35 to ensure that the discretization is entropy 
stable. This is detailed in the following section. 

3.3 Entropy Stable WENO Finite Differences 

It is well known that non-dissipative numerical methods cannot be used to simu- 
late shocks. The primary reason for this is that shocks dissipate energy and non- 
dissipative numerical methods have no mechanism to mimic this. To simulate prob- 
lems with shocks, dissipation needs to be added to the numerical method. There 
are a variety of mechanisms to achieve this, but in this work WENO finite difference 
methods are used. The implementation uses unique formal boundary closures from 
Fisher et al. [10] that satisfy the SBP condition. Stencil biasing mechanics follow two 
papers by Yamaleev and Carpenter [28,29]. The details of the generally applicable 
correction procedure are detailed below. The full implementation details including 
the WENO stencil biasing algorithm throughout the domain are available in Fisher 
et al. [10] and Carpenter et al. [30] for (2-4-2) and (3-6-3) operators, respectively. 1 

The first step to construct a WENO finite difference operator is to cast the 
difference operator in flux form, 

Qf = Af . 

These fluxes that recover the non-dissipative first derivative approximation are called 
target fluxes. The target fluxes are broken into a sum of fluxes on smaller stencils 

1 The nomenclature (2-4-2) signifies that boundaries/interior stencils are second- and fourth- 
order accurate, respectively. The resulting operators are globally third- and fourth-order accurate, 
respectively. 
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of width, p, called candidate stencils , 

n s 

£ = EEE> i = 1,2, - - - , AT - 1, (3.36) 

k= 1 

where n s is the number of candidate stencils needed to describe the target flux, 
f- k are the candidate fluxes, and d k - are the target weights that recover the target 
flux. The candidate stencil width is held constant for all fluxes in the domain. For 
example, fourth-order operators use p = 2, and sixth-order operators use p = 3. 
Figure 2 shows the fourth-order case. The number of candidate stencils needed to 
describe the fluxes fj can vary, as the target fluxes do not all have the same stencil 
size when approaching the boundary. The functional form of the candidate stencils 
depends on the distribution of the flux points, and thus is fully described by the 
norm, V , and the desired order of accuracy. 


Ui—2 1 


Ui 




^i+l 


Ui+2 


Ui + 3 


Si S 2 s 3 s 4 s 5 


Figure 2. The stencil for a WENO scheme with p = 2 and n s = 5 candidate stencils 
is shown. 


WENO works by preventing the interpolated fluxes, fj, from using data across 
discontinuities. This is done by replacing the target weights, cfl) with nonlinear 
weights, 





e 


dj k = d'i 



k = l,...,n s . 


(3.37) 


The functional form of the nonlinear weights relies on the scaling parameter, e, and 
dual stencil-biasing parameters, f and ft. t is a measure of the smoothness over the 
full stencil, 

Tl=t(^0 1 W 2F -') 2 - n T = n,- P . (3.38) 

ft is a measure of the smoothness over each individual candidate stencil, 




p - 1 


£(fe) 2< 


/ dhfjpCj) 

dx l 


2 


(3.39) 


where p k (x) is the unique order (p — 1) polynomial fit of the solution over the 
candidate stencil, Sk- The flux of the WENO scheme is calculated using the formula 

n s 

fj W) = E j = 1, 2, . . . , JV - 1. (3.40) 

k= 1 
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The stencil biasing procedure ensures that stencils containing discontinuities result 
in high relative values of j3 and are assigned negligible weights. The WENO flux 
becomes an interpolation that incorporates only smooth data. Away from disconti- 
nuities, the flux collapses to the target flux. Near discontinuities, the flux transforms 
into an upwind operator. Note that the flux consistency condition on the first and 
last flux point in 2.39 is enforced. 

The WENO method as described has no provable stability properties. We em- 
ploy a limiting procedure to ensure that it satisfies entropy stability. From 3.35, we 
require that 

(»« - »,) T ( fT W) - fT) < o, (3.41) 

where f ( ss ^ ) is the entropy stable WENO flux and is the entropy consistent 
flux. The limiter that guarantees entropy stability is chosen as 

fi ssw) = ir + s (/?> - ip) . * = 

v ' V b 2 + c 2 ( 3 . 42 ) 

b = (w i+1 - Wi) T , c = 1CT 12 , 

where f - " ^ is the WENO flux without the correction described above. The limiting 
process to get the entropy stable property is not unique. It was chosen merely 
because it is smooth with respect to the solution. This same limiting procedure 
works with fluxes other than WENO. It is shown in a later section that this minor 
correction to the WENO operator has a large impact. 

Remark. While this limiting process can be used on other dissipative methods to 
ensure the entropy stability property, we note that the selected dissipative method 
must be cast into telescopic-flux form and satisfy the generalized summation-by- 
parts property. In other words, the dissipative part of the method must have formal 
boundary closures and satisfy the entropy consistency condition. 


3.4 Entropy Stable Viscous Terms 

Using the formalism introduced in Section 2.3.3 combined with the definition of the 
entropy dissipative regularization in Section 2.2, we can define viscous terms such 
that the continuous entropy properties are mimicked by the semi-discrete equation. 
This is accomplished simply by requiring that the discrete viscous fluxes are written 
in dependence on the discrete gradients of the entropy variables, 

0 cw x ) x = T~ l Af (t,) = V 2 (c) w + T P 2 pp, 

V 2 (c)w = V~ 1 (-M(c) + B[c]V)w, ^ 

M(c) = V t V[c\V + K(c), n(c) = Y j N£\c} k N k . 

k = 1 

The accuracy requirements are automatically satisfied. The coefficient matrices [c] 
and [c] are positive semi-definite because they are constructed using block-diagonal 
combinations of positive semi-definite matrices. 
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The contribution of the viscous terms to the semi-discrete entropy decay rate is 


n r 

w T AfM _ Wy gj£]p w _ (X> w ) t T’[c](I>w) — y~](A/'fcw)' r [c] fc (A4w), (3.44) 

k = 1 


where the last term is design-order small, and the last two terms are negative semi- 
definite. Note that no definite properties can be defined if the discrete viscous fluxes 
are constructed based on gradients of primitive or conservative variables. While at 
the continuous level, w x = w q q x , at the discrete level in general we must assume 
Vw ^ w„Du. As in 2.18, only the boundary term will result in growth of the 
entropy, and thus this approximation of the viscous terms is entropy stable. 

This formal definition of high-order entropy stable viscous terms with full bound- 
ary closures appears to be completely new (to our knowledge) . It can be used along 
with the entropy consistent inviscid terms to construct a semi-discretization that 
mimics the properties of the continuous entropy, which is shown in the next section. 


3.5 Entropy Stable Semi-Discretization 

We substitute the entropy stable fluxes for the inviscid 3.42 and viscous 3.43 terms 
into the semi-discrete conservation law 2.43 to get what we define as an entropy 
stable semi-discretization on the finite domain, 

u t + p~iAf (5SW) = + 7 7_1 g 6 . (3.45) 

Note that the boundary conditions imposed with the penalty have not been con- 
sidered. This is because we have found it necessary to specify the conservation law 
and the entropy-entropy flux pair to derive appropriate boundary conditions. In the 
case of the compressible Navier-Stokes equations, the continuous theory for entropy 
stable boundary conditions is lacking, and no satisfactory entropy stable boundary 
conditions have been derived. Instead, linearly stable boundary conditions are used 
in the penalties. The particular boundary conditions used for each conservation law 
are addressed in the next sections. 


4 Applications 

The new methodology developed herein is applied to two sets of conservation laws. 
Developments in Burgers equation are used to illustrate how a fully entropy stable 
scheme is constructed. The methodology is also extended to three dimensions for the 
compressible Euler and Navier-Stokes equations, which are the target application 
for this work. As mentioned above, full entropy stability cannot be achieved because 
of incomplete continuous theory at the boundaries. The entropy stable numerics are 
applied where possible. 
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4.1 Burgers Equation 

Burgers equation is a nonlinear model problem that admits shocks in the inviscid 
limit, 

qt + f{u) x = ef^ v \ /(<?) = y, f {v) = qx, xe[x L ,x R ], te[0,oo), 
q(x L ,t) + \q(x L ,t)\ 

2 9(z.l, t) - eq x (x L ,t) - g L (t) = 0, ^ _ 

q(xR,t) - \q(x R ,t ) | , ,s . x.s n 

2 i) - £q x {XR, t ) + = 0, 

g(x,0) = fifo(x). 


The boundary conditions in 4.1 are constructed such that the entropy corresponding 
to 


(S,F) 


q 2 q 3 

Y’Y 


(vbVO = 


g 2 g 3 

T’ 6 


(4.2) 


only increases with respect to the imposed data and maintains the same form in the 
inviscid limit e — > 0. The entropy variables are w = q. We substitute this definition 
into the entropy decay rate 2.18 and find 


d_ 

dt 



dx 


X L 


eqq x 


XR 


— £ 


J x L 


xr 

/ 


q x q x dx. 


XL 


(4.3) 


Using the boundary condition in 4.1 to replace eq x at the boundaries above, the 
entropy decay rate becomes 


XR 



XL 


X R 

kill 2 k-L I 2 I . f 2 J 

3 q R y q L + qR9R + 9l9l - £ / q x ax, 

XL 


(4.4) 


where qL = q(xL,t) and q R = q(x R ,t). It is immediately clear that when qL = 0 or 
q R = 0, the corresponding boundary terms do not contribute to the entropy decay. 
When qL / 0 or q R ^ 0, the identity, 

~\. { \ 2 -j- 

yz =~2 + Y + Y’ ° >0 ' (4 - 5) 

is used to transform the corresponding positive terms in 4.4, which yields the entropy 
decay rate 


d_ 

dt 



XL 


XL 


1 

2 



+ 



1 v 
v^ 9R ) 

2 if a L 

q R+[Y 





2 , 1 2 , 1 2 
Ql + 2^ 9r+ 2^ 9l - 


(4.6) 
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The coefficients used to complete the square, clr and a R , are arbitrary, and thus it 

is acceptable to set them to values that satisfy the proof, 

2 2 
0<a R <^\q R \, 0 < a L < -\q L \. 

This shows that the only terms in 4.6 that will contribute to growth of the entropy 
are the boundary data terms, ^ 7T[9l an d 2 o~^ 9 r- This result is independent of the 
value of e. 


4.1.1 Entropy Stable Discretization of Burgers Equation 


The inviscid terms in the discretization of Burgers equation are specified according 
to 3.42, 3.9, and 3.10, where the two-point entropy consistent flux simplifies to 


fs{u k , ug ) = - {u k u k + U k Ui + UgUg) 

fa 


(4.7) 


The viscous terms are approximated according to 3.43 with the recognition that 
c = 1. For Burgers equation, entropy stable boundary penalties can be constructed, 
which are equivalent to the penalties derived in [15], 


gb = - ei 


( u i + M / n x \ 

( g ui - e [Du ) 1 - g L J 


+ e N 


UN ~ |«iv| 


u N - e (Du) N + g R 


(4.8) 


The e vectors ensure that only the first and last elements of g & are nonzero. The 
contribution of the penalties to the discrete entropy is 


T 

U gb = - 


+ 


U\ + [till 2 


u( - £U\ (Du) l - uig L 


UN - |«iv| 2 / n \ , 

u N - eu N (Du) n + u N g R 


(4.9) 


Note that the form of the penalties matches the boundary condition specification in 
4.1. The individual discrete elements are combined using 3.45. The semi-discrete 
entropy decay rate is 

^u T Pu < - ^ (u 3 n - it?) + ^ (u% -u\- \u n \u 2 n - M « i ) ^ 10 J 

- eu T Mu + uig L + u N g R . 


Using the relation in 4.5 to transform the positive terms yields the entropy expression 

ji T vs < - eu t mu + T-si + _1_, ’ 

1 


9r 


l - l (Vlo-L)^ - -J=g L 


. OR _ Kjyj 
+ 1 2 3 


* i i a L |^l| \ 2 

U N + ( ~7) g - ) u li 


(4.11) 


2 2 
0<az,<-M, 0<a R <-\u N \- 
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Just as in 4.6, the only terms that lead to entropy growth are an d 2 aZ@R' This 

ensures that the semi-discrete solution will remain bounded as long as appropriate 
boundary data are used. 

4.1.2 Advantages of Entropy Stability 

We show in this section how schemes designed for linear stability may fail at shocks 
and illustrate how entropy stable schemes can overcome such deficiencies. A linear ly- 
stable WENO scheme, termed ESWENO, was developed in Fisher et al. [10] with 
full SBP boundary closures. This scheme is constructed such that when the WENO 
operator is cast in flux form, 

A f (ESW) = (g + 7^ f j q + q t = b , n = n T , (4.12) 

the symmetric dissipation matrix 1Z is positive or negative semi-definite for waves 
with positive or negative propagation speeds, respectively. In comparison with the 
base WENO operator used in this work, the ESWENO operator incorporates an 
additional correction matrix, 

Af (ESW) -Af {W) =K c t, (4.13) 

where 1Z C is a symmetric semi-definite matrix. This correction matrix has been 
shown to be necessary in linear equations where the one-sided stencil biasing me- 
chanics become poorly behaved if no correction is added and result in solutions that 
do not converge. However, it has been observed that the ESWENO scheme exhibits 
overshoots at stationary shocks, which is not surprising because it is based on lin- 
ear theory. This is illustrated by recognizing that the contribution of the energy 
correction term to the nonlinear entropy is 

U T ^ A f(B.sw) _ A f("0) = u T n j = 1 u T ncUu } u = diag(u). (4.14) 

The ESWENO procedure ensures that 1Z C is semi-definite, but the matrix 1Z C U may 
be indefinite and thus may contribute to entropy growth. The disconnect between 
the intended behavior and the actual behavior of the correction term illustrates 
the shortfalls of linear analysis. Instead, if the entropy stable WENO (SSWENO) 
scheme developed herein is applied, the WENO scheme will only dissipate entropy. 
The above behavior is illustrated using a common example defined by 

x L = -1, x R = l, e = 0, t G [0, 2.5], ^ 

9 Lit) = 1, g R {t) = -1, g 0 (x) = -x. 

This problem is discretized using N = 65 points. Figure 3 compares the simulation 
results at t = 2.5 obtained using the WENO, ESWENO, and SSWENO schemes. 
Clearly, the linearly stable ESWENO operator exhibits a dramatic overshoot at the 
shock, while the WENO and SSWENO schemes are monotone. 

While the WENO solution may appear to be the best solution, that is not in 
general what is observed. The SSWENO scheme exhibits better stability properties 


27 



Figure 3. SSWENO, ESWENO, and WENO results for the stationary Burgers shock 
case are compared for N = 65 points at t = 2.5. 


for more complex problems. This is demonstrated with a different test problem with 
time-dependent boundary conditions, 

x L = -l, XR = 1, £ = 0, t E [0, 5], 

1 1 (4.16) 

g L (t) = 1 + - sin(Trf), g R (t) = -1 - - sin(vr t), gi(x ) = -x. 

Multiple shocks meet at x = 0 in this configuration. The problem described by 4.16 
was simulated with N = 65 points and the SSWENO results are compared to the 
WENO and ESWENO results at t = 2.5 (see Figure 4). It is clear that SSWENO 
performs as expected. The overshoots at the shocks disappear, and the one-sided 
stencil biasing mechanics are well controlled. (Figure 4). This problem also demon- 
strates how the WENO scheme without correction degenerates at boundaries. 


4.2 Euler and Navier-Stokes Equations 

The end goal of this work is to simulate the Navier-Stokes equations in a stable and 
efficient manner. Toward this end, the application of entropy stable inviscid and 
viscous terms are detailed in this section. The developments for the Navier-Stokes 
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Figure 4. The solutions using SSWENO, ESWENO, and WENO on N = 65 nodes 
are plotted for the multiple shock problem at t = 2.5. 

equations in this work are limited to the calorically perfect form, 

qt + (f)xi = (f (v)l ) Xi , xen, te[0,oo), 

Bq = gbi x 6 <912, t E [0, oo), (4-17) 

q(x,0) = go(x), x£tt, 

where the Cartesian coordinates, x = (xi,x 2 ,x 3 ) T , and time, t, are independent 
variables and index sums are implied. The three dimensional domain is defined by 
the box 

n = [zf,xf] x [x%,x%] x 

and represents the boundary of the domain. The conservative variables are 

Q = (d, pv \ , pv 2 ,pv- A , pE'f , (4.18) 

where p denotes density, v = (v\,V 2 , v^) T is the velocity vector, and E is the specific 
total energy. The convective fluxes are 

f = ( pvi , pviv i + Sup, pviv 2 + S i2 p, pviV 3 + 5 i3 p, pviH) T , (4.19) 

where p represents pressure, H = E + p/p is the specific total enthalpy, and Sij is 
the Kronecker delta. The viscous flux terms are 

/(")* = (o, Til,Ti2, Ti 3 , TjiVj - qi) T , (4.20) 
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where the shear stress is 


T ij — P Vi)xj + (' Vj) Xi 8ij \ , (4-21) 

and the heat flux is 

q i = -KT Xi . (4.22) 

T denotes the static temperature, and p = p(T), and k = k(T) are the dynamic 
viscosity and thermal conductivity, respectively. The viscous terms can also be 
expressed as 

r r>l (4.23) 

which is a convenient form for entropy analysis. The constitutive relations for a 
perfect gas are 

h = H — -VjVj = c p T, (4.24) 

where c p is the constant specific heat, and 

v=p RT ’ < 4 - 25 > 

where R u is the universal gas constant and MW is the molecular weight of the gas. 

The speed of sound for a perfect gas is 

c = yfyKT, 7 = ~^~n' ( 426 ) 

Cp it 

In the entropy analysis that follows, the definition of the thermodynamic entropy is 
the explicit form, 

<427) 


4.2.1 Entropy Analysis 

In the Navier-Stokes equations, entropy stability is not the same as full nonlinear 
stability. However, entropy stability gives a stronger stability estimate than linear 
energy stability, and in many ways is easier to apply. In this section, the continuous 
entropy stability is conducted first to illustrate the entropy characteristics of the 
governing equations at the continuous level. Then, spatial operators are derived 
that enable these continuous properties to be mimicked, which is shown through the 
semi-discrete entropy analysis. 

The entropy-entropy flux pair for the Navier-Stokes equations is 

S = — ps, F l = —pviS, (4.28) 

and the potential-potential flux pair is 

tp = pR, il> 1 = pviR. (4.29) 
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Note again that the mathematical entropy has the opposite sign of the thermo- 
dynamic entropy. To avoid confusion, herein entropy refers to the mathematical 
entropy unless otherwise noted. The entropy variables using the pair in 4.28 are 


w = 


Si 



v i v j 

2 T 


Vi V 2 

rj~i ? rj~i ' 


V3 

T' 



(4.30) 


and can be shown to have a one-to-one mapping with the conservative variables as 
long as p,T > 0. Expressly: 

( T S qq ( T > 0, VC + 0, p,T> 0. 


This restriction is what makes the entropy proof fail to be a true measure of nonlinear 
stability. Another mechanism must be employed to bound p and T away from zero to 
ensure stability. This is not considered in the present work, but interested readers 
can refer to the work by Zhang and Shu [31]. The full derivation of the entropy 
variables and symmetrizing matrices are detailed in Appendix B.l. 

The entropy equation is found by premultiplying the Navier-Stokes equations 
with the transpose of the entropy variables, 

S g qt + S q (f) Xi = S t + F l x . = w T ( Cijq w w Xj ) Xi ^ 

= (w T Cijq w w Xj ) xi - w x .Cijq w w Xj , 

where the viscous terms satisfy 

C-ij^w = &ij = Cji, Cxi^ijCxj Ci 0, VC, (4.32) 


requiring 


T > 0, p(T) > 0, k(T) > 0. 


This requirement and the symmetric coefficient matrices, tjj . are derived in Ap- 
pendix B.2. The total entropy decay rate is found by integrating 4.31 over space, 


^ J SdV= J (- w T c ijWx . - F) OS 1 - J w T Xi c ijWx . dV. 

n an n 


(4.33) 


Dutt [32] used the above equation to derive well-posed boundary conditions. His 
procedure has not gained wide acceptance [33]. In this work, boundary conditions 
are still imposed based on linear stability theory. 


4.2.2 Discretization Notes 

To facilitate the extension of the entropy stable methods to the three-dimensional 
equations, we define the three-dimensional nomenclature and examine the general 
form of the semi-discretization. The semi-discrete form of 4.17 is 

U * + E ( r - J{v)l ) = E (gb + g>) • ( 4 - 34 ) 

i= 1 i=l 
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The solution vector is ordered as 


u = (u(* ( i )( i )( i)fXw )(2) ) T , ' ' 


u (®(Afi)(Ar 2 )(Ar 3 )) T ) 


T 


The roman superscript indices on the flux in 2.43 indicate the direction of the flux, 
and parenthetic superscripts indicate the type of flux, i.e. V for viscous, W for 
WENO, S for entropy consistent. The multi-dimensional operator nomenclature is 
defined as follows: 


V Xl = (Vni ® In 2 ® In 3 <8> J5) , T ) x 2 = {In-l ® Vn 2 ® In 3 ® h) , 
V X3 = ® In 2 <8> Vn 3 ® h) , T’ Xl a; 2 = ("Pjvi <8 Vn 2 <8 In 3 8> is) , 

V x , 3,3 = (TVi <8 iw 2 <8 TV 3 <81/5), 7>x 2a * = (iiVi <8 Vn 2 <8 T’wg <8 1.5 ) 


(4.35) 


P = P 


X 1 X 2 X 3 


= (Vn! <8 Vn 2 <8 iV 3 (8> Is) 


where <8> represents the Kronecker product; Vn represents the one-dimensional norm 
defined with N solution points; In is the one-dimensional identity operator defined 
with N solution points; and the other operators, £>, Q, A,A4,7Z, extend in an iden- 
tical fashion. Note that transposition and inverse operations applied to these opera- 
tors are distributive. When applying these operators to the scalar entropy equation 
in space, a hat will be used to differentiate the scalar operator from the full vector 
operator. For example, 


V = (TVg <8 Vn 2 <8 Vn 3 ) ■ 


(4.36) 


Expressly, 4.35 states that one subscript on an operator means the operator applies 
in one direction, two subscripts apply in two-directions, and three subscripts or 
boldface apply in all three directions. This nomenclature is used throughout the 
following analyses and definitions. 


4.2.3 Entropy Stable Spatial Discretization 


The inviscid terms in the discretization of the Navier-Stokes equations are calculated 
according to 3.42, 3.9, and 3.10, where the two-point entropy consistent flux was 
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derived by Ismail and Roe [8] 


f J s (Ui, U i+ 1) = (pVj, f>VjVl + Sjip , pVjV 2 + Sj2Pi pvjvs + 5 j3 p, pvjH 
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v = 


Vi + 1 Pi i Pi+l 

+T y/Ti+l A +T +7i+i 

1 . I ’ P ~ 1 . 1 ’ 


log 


h = R- 


VTjpi 


y/Ti 


i+lPi+l 
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V+ + +7+7 


VR +1+1 
VTiPi + \/7i+lPi+l 


V (v^i + W^iPi VTi+ipi+i) (4.37) 
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7+1 


log ( \/% i 


H = h+ -vpvi, p = 


7 - 1 




log 


+ 


Pi 


^~i+i Pi+l / \ y/Ti ■yJ'T'i - (_i y y 

1 ' (-v/TiPi - v/Tj+i/jj+i) 


2 (log(+Rpi) - log(Vr i+ i/9 i+ i)) 


This somewhat complicated explicit form is the first entropy consistent flux for 
the convective terms with low enough computational cost to be implemented in 
a practical simulation code. Previously, Tadmor [5] derived an entropy consistent 
flux form that required integration through phase space, but this was deemed too 
expensive to be practical. 

The resulting entropy stable WENO scheme for the Navier-Stokes equations 
satisfies 


w 


T *J SSW * > 4 d 


(4.38) 


in each direction. This new version of WENO is more robust than ESWENO for 
the Navier-Stokes equations and easier to implement. The current major limitation 
is that it has only been derived for the calorically perfect equations. 

The definitions in 4.35 are used to extend the viscous operator in 3.43 to a 
system of equations in three dimensions. Recall that to demonstrate the entropy 
decay resulting from the viscous terms, the discrete viscous fluxes must be written 
in terms of entropy variable gradients, 


( CiiW Xi )xi — ^XiXi {pii ) W + Tp2ppi 

'DxiXii&ii ) w = 'Pxi (—■M-x i (cu) l3 Xi [cii\D Xi ) W, 

n r 

■Mxi(cii) = 'Dxi'Pxi \c-ii\Dxi + 77^ (cjj) , 'R-x i (cii) = ^ ^ ^kxi [<+] k^kxj ■ 

k= 1 


(4.39) 


The definition above only considers the approximation of the diagonal viscous terms, 
where two derivative operators are applied in the same direction. Cross terms can 
be incorporated into the same form using the SBP property, Q = B — Q 1 , yielding 
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the full viscous flux gradient, 


(c-ijW X i)xj — 7e? 

VxiXjiCij) w = P,. : l {-'D^.T> x .[c i j]V x . - dijR.r: {('ij ) + B Xi [c ij ]V Xj ) w, 

n r 

'R'Xj(Cij) = -^kxH [Cij\k/^ kx i • 
k = 1 


(4.40) 


We recast the viscous gradient in 4.40 into a flux difference form. The equivalent 
viscous flux differences are defined by 


3 n r 

Aijf = — w — yy A/fo!. [ciij^A/^x^w. (4.41) 

€=i £=i 

— (l?) 

The derivation of the actual flux form, f , is somewhat complicated and long, 
including two parts. The procedure to calculate the contribution of the diagonal 

— ('ll) . 

viscous terms to f is detailed in Appendix A. 2.1. The procedure to calculate 

— (.U ) 

the contribution of the cross terms to f follows from Appendix A. 1.1, where the 
transverse gradients, calculated using SBP finite differences, are simply interpolated 
to the flux point using the high-order interpolation formulas. 

We calculate the contribution of the viscous terms to the semi-discrete entropy 
decay rate by premultiplying the flux difference in 4.41 by 'w T 'P, 


w r v yy ■ = yy yy wi , pv Xl i B Xi [c u }v Xf w 


i= 1 


i= 1 i= 1 
3 3 


yy yy w ) T @>x t w) 

i = 1 1=1 
n r fir 

yy yy (a/'&.w ) 1 w~^[cn] t {M tXi w) 


(4.42) 


i = i e=i 


Reorganizing the gradient terms into a long vector, we show in Appendix B.2 that 
the latter two terms in 4.42 are always negative semidefinite as long as temperature, 
viscosity, and thermal conductivity are positive everywhere in the domain. 

The multidimensional form of the semi-discrete entropy decay is found by com- 
bining the contributions from the entropy stable WENO flux in 3.42 and the entropy 
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stable viscous terms in 4.42. The resulting semi-discrete entropy decay rate is 



1=1 


1=1 


3 3 


+ ™ T ' jyP Xi B x i [Ci£Px t W 


i = i e = i 
3 3 


(4.43) 


- X! 2 ( p *i w ) T (^ w ) 


i=l £=1 


I t/'p I 

- 5^ X) (A/'& i w) T 7 7 'P“ 1 [cjj]^ (JV&.W) . 


i=l £=1 


The last term is design-order small and vanishes as the grid is refined. The other 
terms are analogous to the terms in 4.33. For high-order finite-difference (HOFD) 
methods, this is the first known implementation that satisfies this property. 


4.2.4 Energy Stable Boundary Conditions 

All the problems we have studied in this work have been limited to open boundary 
conditions, which are implemented in the form of Svard et al. [22] Extensions of these 
energy stable boundary conditions to walls are detailed in Berg and Nordstrom [24]. 


5 Accuracy Validation and Robustness 

The accuracy and robustness of the algorithms developed herein are tested using two 
smooth and two discontinuous problems. The smooth problems are the propagation 
of an isentropic vortex, and the propagation of the viscous shock. Both problems 
demonstrate the design-order convergence of the new entropy consistent formulation 
and of the SSWENO scheme. The final two problems (the Sod and Lax shock tube 
problems) have discontinuous solutions, and test the efficacy of the entropy stable 
correction terms used in the SSWENO formulation. It is demonstrated that the 
corrections do not adversely affect the desirable stencil biasing properties of the 
baseline WENO scheme. 


5.1 Isentropic Vortex 

The isentropic vortex is an exact solution to the Euler equations that has proven to 
be an excellent test of the accuracy and functionality of the inviscid components of 
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a Navier-Stokes solver. It is fully described by 

f{x, y, z, t) = 1 - {x-xq- Uoc cos (a) tj 2 + {y - y 0 - sin(a) t ) 2 
T(x,y,z,t) = 1 - elM^^exp(f(x,y,z,t)) , p(x,y,z,t) = 

u(x,y,z,t ) = I/ooCOsO) - e v y ~ Vo ~ U ^ sm ( Q k CX p 

v(x,y,z,t ) = I/oo sin(a) - e v x ~ x °~ U ^ cos(Q)f exp 
w(x,y,z,t) = 0. 

In this study the values C/qo = e v = 5.0, Moo = 0.5, and 7 = 1.4 are used. 


(5.1) 


5.1.1 Optimal Accuracy: A Periodic Cartesian Grid Test Case 

Theorem 3.1 proves that design order entropy consistent fluxes can be constructed 
using a linear combination of two-point entropy fluxes. A critical assumption used 
in the proof is that the two-point non-dissipative fluxes satisfy Tadmor’s integral re- 
lation given in equation 3.10. Herein, the non-dissipative Euler fluxes of Ismail and 
Roe [ 8 ] are used, which to our knowledge have not been shown to satisfy Tadmor’s 
integral relation. The first study is designed to detect the presence of order reduc- 
tion resulting from the use of the Ismail and Roe fluxes in the entropy consistent 
formulation. A periodic domain is used to eliminate the order reduction associated 
with p th -order boundary closures in the diagonal norm operators, thus enabling the 
potential of fully 2 p i?t -order convergence. 

The Cartesian grid test case is described by 

x € (—15, 15), y € (—15, 15), (x 0 ,y 0 ) = (0,0), a = 0.0, t> 0. 

Although periodic boundary conditions were used, note that the actual propagating 
vortex solution is not periodic. For this domain and the chosen vortex strength, the 
deviation from periodicity is exponentially small. Discrete error is compared with 
the exact solution after one full flow through time. 

Three grid resolutions are examined using both the fourth- and sixth-order en- 
tropy consistent target finite difference schemes and the entropy stable WENO 
schemes. The error decay, shown in Tables 1 and 2, asymptotes towards the de- 
signed rate in each case (i.e. , fourth-order and sixth-order respectively). The error 
for the SSWENO scheme approaches the target central scheme as the grid is refined. 
The increase in the SSWENO error constant is more pronounced in the sixth-order 
case, although still within about a factor two of the target entropy consistent FD 
scheme. 

The study provides evidence that the the non-dissipative Euler fluxes of Ismail 
and Roe [ 8 ] do not degrade the formal accuracy of the high-order entropy consistent 
fluxes. 


5.1.2 Finite Domain Cartesian Grid Test 

Theorems 3.1 and 3.2 provide a recipe to construct design-order entropy consistent 
fluxes on a finite domain, including points with biased stencils close to each bound- 
ary. Having established fully 2p th -order convergence for the interior operator, now 
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Table 1. Error convergence is shown for the fourth-order simulation of the periodic 


isentropic vortex. 



Entropy Consistent FD 

Entropy Stable WENO FD 

Resolution 

L 2 error 

L z rate 

L°° error L°° rate 

L z error L z rate 

L°° error L°° rate 

60 x 60 

1.05 

- 

0.38 

0.88 

0.29 

120 x 120 

7.60e-02 

3.79 

3.19e-02 3.58 

8.85e-02 3.32 

3.38e-02 3.12 

240 x 240 

5.01e-03 

3.92 

2.11e-03 3.92 

5.97e-03 3.89 

2.81e-03 3.59 


Table 2. Error convergence is shown for the sixth-order simulation of the periodic 


isentropic vortex. 



Entropy Consistent FD 

Entropy Stable WENO FD 

Resolution 

L 2 error 

L 2 rate 

L°° error L°° rate 

L 2 error L 2 rate 

L°° error L°° rate 

60 x 60 

0.31 

- 

0.13 

0.30 

0.11 

120 x 120 

6.05e-03 

5.67 

2.50e-03 5.69 

1.20e-02 4.63 

5.58e-03 4.34 

240 x 240 

9.76e-05 

5.95 

3.81e-05 6.04 

2.46e-04 5.60 

1.83e-04 4.93 


a finite-domain Euler vortex test case is used to establish convergence rates for the 
boundary closures when using a diagonal norm, finite-domain operator. 

The finite domain Cartesian grid test is described by 

x <E (-5, 5), y G (-5, 5), (x 0 , yo ) = (0, 0), a = 0.0, t > 0. 

Three different grid resolutions are examined, with the vortex located precisely on 
the boundary when the error measure is evaluated. This measures the effect of the 
boundary closure, the penalty boundary condition, and the interior scheme. Both 
the (2-4-2) and (3-6-3) entropy stable WENO finite difference schemes are evaluated. 
The designed order of accuracy (p + 1) is observed for all algorithms as indicated 
in Tables 3 and 4, in accordance with the classical theorem of Gustafsson [34], The 
reduction in order from 2p to p + 1 relative to the periodic test case results from the 
order p boundary closures. No explanation is currently available for why the (3-6-3) 
scheme appears to converge at a super-optimal rate 


Table 3. Error convergence is shown for the (2-4-2) simulation of the finite domain 
isentropic vortex. 



Entropy Consistent FD 

Entropy Stable WENO FD 

Resolution 

L z error 

L z rate 

L°° error 

L°° rate 

L z error 

L z rate 

L°° error L°° rate 

32 x 32 

5.37e-02 

- 

3.17e-02 

- 

5.44e-02 

- 

3.49e-02 

64 x 64 

6.14e-03 

3.13 

3.51e-03 

3.18 

6.26e-03 

3.12 

3.51e-03 3.31 

128 x 128 

7.43e-04 

3.05 

4.63e-04 

2.92 

7.46e-04 

3.07 

4.72e-04 2.90 

256 x 256 

9.01e-05 

3.04 

5.95e-05 

2.96 

8.39e-05 

3.15 

5.52e-05 3.10 


5.2 Viscous Shock 

The entropy consistent Navier-Stokes formulation developed herein, relies on newly 
derived narrow stencil viscous terms. In principle, design order convergence should 
be achieved for all terms including the nonlinear viscous terms in the energy equa- 
tion. The convection of the viscous shock (an exact solution to the Navier-Stokes 
equations) is used to test this assertion as well as the functionality of the viscous 
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Table 4. Error convergence is shown for the (3-6-3) simulation of the finite domain 


isentropic vortex. 



Entropy Consistent FD 

Entropy Stable WENO FD 

Resolution 

L z error 

L 1 rate 

L°° error 

L°° rate 

L z error 

L z rate 

L°° error L 

00 rate 

32 x 32 

7.33e-02 

- 

3.65e-02 

- 

7.32e-02 

- 

3.73e-02 

- 

64 x 64 

4.30e-03 

4.09 

2.10e-03 

4.12 

4.95e-03 

3.89 

3.03e-03 

3.62 

128 x 128 

1.94e-04 

4.47 

2.08e-04 

3.34 

2.53e-04 

4.29 

3.80e-04 

3.00 

256 x 256 

7.72e-06 

4.65 

1.79e-05 

3.54 

1.09e-05 

4.54 

3.09e-05 

3.62 

512 x 512 

2.83e-07 

4.77 

1.01e-06 

4.15 

4.19e-07 

4.70 

1.71e-06 

4.18 


terms in the solver. The nonlinear convection and viscous terms are perfectly bal- 
anced in this computation, and thus the shock thickness remains constant, and 
the shock front is simply advected. The derivation of this problem is available in 
Fisher [35]. 

The convergence rate for the viscous shock is evaluated on a Cartesian grid, 
described by 

x € (-1,1) x (-0.5, 0.5). 

The shock flow is rotated a = 20° with respect to the x-axis. The shock profile 
is initially located at x s = —0.5 with respect to the origin and is simulated until 
t = 0.25. The Reynolds number was Re = 10, and the reference Mach number was 
M = 2.5. Table 5 shows the errors for the Cartesian grid. Because the Navier- 
Stokes equations are incompletely parabolic in type, the design order convergence 
rate for the (2-4-2) and (3-6-3) operators should be 3 and 4, respectively. Once 
again, super-optimal convergence is observed for this test case. No explanation is 
currently available for why the schemes converge at a super-optimal rate. 

All smooth test cases demonstrate the validity of the Theorems 3.1 and 3.2. 
Furthermore, the SSWENO formulation achieves design order in all cases with only 
marginal loss in accuracy relative to the baseline target operator. 


Table 5. Error convergence is shown for (2-4-2) and (3-6-3) simulations of the viscous 
shock. 



Entropy Stable(2-4-2) 

Entropy Stable (3-6-3) 

Resolution 

L z error 

L z rate 

L°° error 

L°° rate 

L z error 

L z rate 

L°° error 

L°° rate 

48 x 24 

1.74e-03 

- 

5.60e-03 

- 

3.36e-04 

- 

1.57e-03 

- 

96 x 48 

1.10e-04 

3.97 

4.24e-04 

3.72 

7.62e-06 

5.46 

6.51e-05 

4.59 

192 x 96 

7.30e-06 

3.92 

3.38e-05 

3.65 

2.53e-07 

4.91 

3.42e-06 

4.25 


5.3 Shock Tube Problems 

All WENO schemes including those used herein [10,28,29], incorporate elaborate 
stencil biasing mechanics to achieve nearly monotone solutions in the vicinity of 
shocks and other discontinuities. It is critically important that the additional dis- 
sipation provided by the entropy stabilization terms does not contaminate the de- 
sirable attributes of the baseline WENO scheme. To date, all evidence suggests 
that this it indeed the case (including test cases run on non-Cartesian meshes on 
curvilinear grids, but not reported here). The shock tube problems shown below 
are representative of the observed behavior of the SSWENO scheme. 
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(a) Sod Shock Tube, t = 0.2 (b) Lax Shock Tube, t = 1.3 

Figure 5. Shock tube solutions are plotted for the entropy stable WENO methods 
developed in this work and compared to reference solutions. 

5.3.1 Sod Shock Tube 

Sod’s shock tube problem evaluates the behavior of a numerical method when a 
shock, expansion, and contact discontinuity are present. Of particular interest is 
additional observed smearing in the shock and contact, or oscillations at any of 
the discontinuities. These would be indications of adverse side-effects of the entropy 
stabilization term, since the baseline operator performs adequately on this problems. 
The domain is 

i 6 ( 0 , 1 ), y G (— 0 . 1 , 0 . 1 ), t > 0 , 

and is initialized with 

, \ f 1 x < °- 5 > ( \ / 7 

= j 1/8 *> 0 . 5 , p(w) = ( 7/10 

u(x,y,z) = 0 , v(x,y,z) = 0 , w(x,y,z) = 

where 7 = 7/5. The problem is simulated with the (2-4-2) and (3-6-3) entropy stable 
WENO operators using N = 100 uniform cells. The solution is plotted for t = 0.2 
in Figure 5(a). The solutions do not exhibit oscillations and the shock smearing is 
nearly equivalent between the two schemes, with slightly less diffusion observed in 
the (3-6-3) scheme. All profiles are essentially equivalent to those obtained with the 
baseline WENO operator. 

5.3.2 Lax Shock Tube 

Lax’s shock tube problem is used to show that no entropy problems are observed 
using the current methodology and that the correct shock location is observed. The 
reference solution uses N = 800 points with the (2-4-2) WENO operator. 


x < 0.5, 

x > 0.5, ( 5 . 2 ) 

0 , 
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The simulated domain is 



x G ( 

-5,5), y 

G (-0.5, 0.5), 

t > 0, 


initialized with 






p{x, V,z) = | 

0.445 

x < 0.0, 

P(x,y,z ) = • 

f 3.5287 

x < 0.0, 

0.5 

x > 0.0, 

[ 0.5717 

x > 0.0, 


f 0.698 

x < 0.0, 



O 

II 

u(x, y, z) = 

\ 0.0 

x > 0.0, 

v (x, y, z) = 

0, w(x, 


where again 7 = 7/5. The simulation used N = 200 uniform cells and the solution 
is plotted in Figure 5(b) for t = 1.3. Again the solutions do not exhibit oscillations. 

6 Conclusions 

High-order entropy stable finite difference methods have been derived for conserva- 
tion laws on finite domains. These methods use formal boundary closures to satisfy 
a generalized summation-by-parts property. A new method was developed to ensure 
that dissipative inviscid approximations for shock capturing are guaranteed to be 
entropy stable. This was shown for WENO but is more generally applicable. A new 
high-order entropy stable viscous approximation method was also developed using 
narrow-stencil finite difference operators. These new methods were shown to fix 
deficiencies observed in linearly stable methods and do not adversely affect shock 
capturing capability in Burgers equation or the Euler equations. The approach was 
shown to preserve design-order accuracy of diagonal SBP finite difference operators 
on simulations of the Euler equations and the Navier-Stokes equations. 
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Appendix A 


Summation- by-parts operators— (2-4-2) 


The (2-4-2) family of finite difference operators used in this work is specified in 
this section. This operator set enables the Navier-Stokes equations to be simulated 
with third-order accuracy. Another important note is that energy and entropy 
stability proofs require a common diagonal V norm to be used for all terms in 
an approximation. The number of special boundary points for this operator is 
s = 2p = 4. With this specification, the diagonal V norm for the (2-4-2) operator is 
unique, 


V 


diag 


17 59 43 49 49 43 59 17' 

48’ 48’ 48’ 48’ 1 ’ 48’ 48’ 48’ 48 * ‘ T ’ 


(Al) 


where V is an ( N x N ) matrix, and thus has N elements on the diagonal. This form 
follows immediately from the accuracy requirement and form of the first derivative 
operator. 


A.l First Derivative 


Recall that the summation-by-parts operator approximating the first derivative has 
the form, 

V = V~ 1 Q. 


When using a diagonal norm, V, and s = 2p boundary points, this operator is 
unique. The matrix elements at the left boundary have the structure, 


/ _1 59 _j_ _j_ 

/ 2 96 12 32 


Q = 


i 

l 6 

¥ 

32 

0 

0 


59 

96 

1 

12 

1 

32 

0 

0 

0 

0 

59 

96 

0 

0 

0 

0 

59 

96 

0 

0 

59 

96 

59 

96 

0 

1 

3 

0 

1 

12 

0 

0 

0 

1 

2 

0 

2 

1 

12 

3 

3 

l: 2 

0 

0 

1 

12 

2 

3 

0 

2 

3 


0 

0 

0 

0 

0 

_ J_ 

12 


(A2) 


Q is an (N x N) matrix and elements corresponding to the right boundary have the 
property, 

Q(i,j) = ~Q(N+i-i,N+i-j), 1 <i,j<N. (A3) 

Each interior row has the form 


1 2 _ 2 1 

Q(i,i— 2) "pi ’ ^(*,1— 1) 1) Y’ 9(i,i+ 2) ’ 

i = s + 1, . . . , N — s. 


(A4) 
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A. 1.1 Flux Form 


The first derivative SBP operator above can be recast into the flux form, 

Qf = Af = AXf . 


(A5) 


The interpolation matrix, X, is an (N + 1 x N ) matrix where each row sums to 1. 
The matrix elements near the left boundary are 

\ 


X = 


( 1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

2 

59 

1 

12 

1 

32 

0 

0 

0 

0 

0 

11 

f 

32 

96 

0 

17 

f? 

32 

1 

T 52 

12 

0 

1 

12 

0 

0 
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0 

0 

0 

0 

0 

1 

12 

7 

12 

7 

12 

1 

12 

0 

0 

0 

0 

0 

0 

1 

12 

7 

12 

7 

12 

1 

12 

0 

0 

0 

0 

0 

0 

1 

12 

7 

12 

7 

12 

1 

12 

0 


V : 


(A6) 




The flux consistency condition is imposed in the first row of X. The right boundary 
terms satisfy the matrix property, 


h (ij) = h (N-i,N+i-j), 0<i<N , 1 <j<N, 


(A7) 


where h(jj) is the element in row i and column j in X. Note that for the interpolation 
matrix the first row is indexed at 0 to be consistent with the flux point nomenclature 
used throughout this document. The interior elements have the formula, 


T __1 T 7 _l_ T __1 

1) 22’ ^(v) "-(1,1+1) 1 r) > ’ l (i,i+2) 


12 ’ 
i = s, , 


12 ’ 


12 ’ 


(A8) 


, n — s. 


A. 2 Variable Coefficient Second Derivative 


To find the matrices specified in 2.37, A fk and [i?]*,, only i9(x) = 1 must be consid- 
ered. Then, for variable coefficients and nonlinear problems, the structure of \&\k 
ensures the accuracy and stability of the full operator. The structure of Ad(l) is a 
pentadiagonal matrix with (s x s) block matrices superimposed on the diagonal at 
the two boundaries. The first constraint is to require At(l) to satisfy the optimal 
accuracy, 

<Mx) = V- 1 (~M( 1) + BV) cj) + 75-4-2. 

At fourth-order, this leaves one free parameter in the matrix Al(l). The remainder 
is constructed using the matrices: 


M = 


-1 

3 

1 

CO 

1 

0 

0 

0 

0 


0 

0 

-1 

3 

-3 

1 

0 

0 

0 


0 

• o 

0 

-1 

3 

—3 

1 

0 

0 


0 

• • o 


0 

0 

-1 

3 

—3 

1 

0 

0 

0 


0 

0 

0 

-1 

3 

-3 

1 

0 

0 


0 

0 

0 

0 

-1 

3 

-3 

1 


(A9) 
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which is a (N — 3 x N ) matrix, 






( 1 -4 

6 

-4 

1 

0 

0 


0 

0 


0 ^ 







0 1 

-4 

6 

-4 

1 

0 


0 

0 


0 







0 0 

1 

-4 

6 

-4 

1 


0 

0 


0 





A /2 

= 












5 


(A10) 





0 ... 

0 

0 

1 

-4 

6 


-4 

1 

0 

0 







0 ... 

0 

0 

0 

1 

- 4 


6 

-4 

1 

0 







\ 0 ... 

0 

0 

0 

0 

1 


-4 

6 

-4 

1 



which is a 

(N- 

-4 

x N) matrix, 

and 












/ 

-1 

5 

-10 

10 

-5 

1 

0 


0 


0 

0 


0 

\ 




0 

-1 5 

-10 

10 

-5 

1 


0 


0 
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0 





0 

0 

-1 

5 

-10 

10 

-5 


1 


0 

0 


0 



A3 — 

















, (All) 
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0 

0 

-1 

5 

-10 


10 


-5 

1 

0 

0 
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0 

0 

0 

-1 

5 

- 

-10 

10 

-5 

1 

0 
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0 

0 

0 

0 

-1 


5 


-10 

10 

-5 

1 

) 


which is 

a 

(N 

- 5 

1 x N) matrix 

. If only 

these 

matrices are 

used to 

construct the 


remainder, a proof that the remainder is positive semi-definite will not be possible. 
We find that an additional ( N x N) matrix, A/4 is required with only the first and 
last rows containing nonzero terms, 


= — ai + a 2 - a 3 , ra|^ 2 ) = 3a i _ 4 a 2 + 5 a 3 , 

= — 3 ai + 6a 2 - 10 a 3 , = ai — 4 a 2 + 10 a 3 , 


(A 12 ) 


~ ( 4 ) - 

n (l, 5 ) = a 2 - bO- 3 , 


~ ( 4 ) _ 

n (l,6) a3j 


n {N,j) 3 I; 2 , , N. 


This introduces nonlinearity into the derivation of the operator, which is addressed 
differently in our methodology below. The form of the variable coefficient matrices 
is 

1 )(i,i) = i&i+i + $1+2 ) , * = 1 , 2 , . . ., N - 3 , 


(fcM = -c^Vi+ 2 , i = 1 , 2 ,..., N - 4 , 

3 = ~Pp\ (#i + 2 + ^i+ 3 ) , i = 1 , 2 , . . . , N — 5 , 
([^ 4 )m = ^> i = 1 , 2 ,..., AT, 


(A 13 ) 


and for the operator to be provably stable, all coefficients, cp < 0 . For $(x) = 1 , 
only these coefficients are left in [d] ^ . The coefficients, cp can be solved for linearly 
as nonlinear functions of the free parameter left in A 4 (l), ai, a 2 , and a 3 , such that 
the remainder term recovers the narrow form of A 4 (l) with the accuracy constraints 
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satisfied. The resulting coefficients, organized into vectors are: 


c 


(i) 


/ 


if - aia 2 + aia 3 + |to( 2) 4) - 
aia 2 - 2aia 3 - 
a ! a 3 “ 3528 

I¥ 

Jt 

18 


181507 

1719312 


\ 


V « ‘f 


18 185 

° ia 3 3528 


aia 2 - 2aia 3 - 4++ 
aia 2 + aia 3 + gTO( 2)4 ) - 


181507 

1719312 


J 


(A14) 


and 


/ — dia 2 + o 2 — o 2 a 3 — 


c (2) = 


30568 

-(2tti - 0 2 )o 3 - J2 


1^4 

I|4 

I|4 

lj4 

T|4 

l44 


1 

' 144 
1 

' 144 
1 

' 144 
1 

' 144 
1 

' 144 


(2o 3 0 2 )o 3 2g2 


\ — aia 2 + a 2 — a 2 a 3 — 


50568 


C< 3 > 


/ (ai - 02)03 +01+352 \ 
0 
0 
0 
0 
0 
0 


0 

0 

0 

0 

0 

0 

V (ai - a 2 )a 3 + a§ + 355 / 


(A15) 


To eliminate the remaining nonlinear terms such that c>p < 0, a nonlinear 
optimization problem is solved, where the functional is 

c^) +5 + cf\ (A16) 

i j 

(j) 

which is very small as long as all c\ are negative or zero. The functional is mini- 
mized using a conjugate gradient method in the open-source computational rnathe- 
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matics package SAGE [36] 


16815244 _ _ 88998127 

m(2 ’ 4) ~~ 410099621’ 0l ~ _ 304807400’ 
6823462 __ _ 20751280 

m ~ 373821039’ “ 3 “ ~ 551691433' 

The constant coefficient narrow stencil matrix has the form, 


M( 1) 


/ -§"7(2,4) + f "1(2,4) -fl - m (2,4) + A 5>"(2,4) + ® 0 

m (2,4) - || -3"»(2, 4) + |§ 3r71 (2 4) - §| “"7(2,4) 0 

-»"(2,4) + T2 3 111(2, 4 ) - I 3"l(2,4) + f§ 171(2,4) “if h 

l m (2,4) + ® -"1(2,4) "1(2,4) -§ “§"7(2,4) + § 

0 0 A -4 5 

u u 12 3 2 

V 

m ^(JV+1— i,N+l— j) j 1)2 ,... , A r . 


(A17) 
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While thus far only the constant coefficient case has been considered, the variable 
coefficient viscous gradient operator has been fully specified. It is noted, however, 
that the number of special boundary points has increased from s = 4 to s$ = 6 
because the form used here is not optimal for the boundaries. This was deemed 
appropriate for the time being but will be revisited in the future. The accuracy of the 
variable coefficient gradient is guaranteed by using only undivided third derivative 
approximations and above. The lowest-order remainder term is O(Sx^). Higher- 
order methods can follow in a similar way. 


A. 2.1 Flux Form 

The operators above are inconvenient for implementation. Instead, the viscous flux 
gradient is calculated using the flux form, 

V 2 (#)v = V~ 1 Af {v \ (A18) 

where f is consistent with the viscous flux to design order. To calculate the fluxes, 
a coefficient array is needed for each flux point, and a three-dimensional array is 
required to describe the full operator. To reduce confusion, a compressed sparse 
row structure is used below to show coefficients. The pointer for the beginning of 
each row is denoted by 2 . The two column pointers corresponding to the coefficient 
and the variable are denoted by o^ 1 -* and o^ 2 \ respectively. Note that there are two 
because of the three-dimensional coefficient array. The coefficients are denoted by 
b. The form of the viscous flux is 

N N z i+ i-l 

fi = X/ XI C(i,j,k)'&j v k = be$ 0 (.i)V o (2), 0 <i<N. (A19) 

j= 1 k = 1 i=Zi 

Note that the structure of the coefficients is such that 

= C(N+i-i,N+i-j,N+i-k)> 0 <i<N, 1 < j, k < N. (A20) 
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The left boundary and first interior points are described by 


z = (1, 5, 27, 49, 75, 100, 130, 144) , 


(A21) 


and 
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Appendix B 


Navier-Stokes Equations— Supplemental Details 


B.l Derivation of Entropy Variables 

Recall that the entropy variables are defined as w T = S q . While this can be calcu- 
lated directly for the calorically perfect equations, in general it is more convenient to 
use the primitive variables to aid in the derivation. To accomplish this, the standard 
differential relations for pressure, enthalpy, internal energy, and entropy are used, 


dp = pRdT + RT dp, 


dh = 


7 R 

7-1' 


dT , de = 


R 


7- 1 


-dT 


RdT 


dh dp 

S = T ~ w = (7 - 1)T 


-jA 

p 


A set of primitive variables, 

v = (p,v 1 ,v 2 ,v 3 ,T ) T , 

is selected for the derivation. The expansion S q = S v v q is used, where 

pR 


S v = [R-s, 0,0,0,- 


is combined with 

/ 


m Vl 
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Rp 


0 

0 

0 

1 

p 

(7~l)v 3 

Rp 


0 \ 
0 
0 
0 

1=1 / 
Rp ) 


yielding 


c _ 1 _ I p _ _ v kVk Vl V2 V3 _ ^ 

^ l ^ yp 5 p 5 p 5 p 5 p 


A similar procedure is followed to find the matrix w q = S qq = w v v q , where 


w v = 
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p 

0 
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(B2) 


(B3) 


(B4) 


(B5) 
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To illustrate why p, T > 0 are required for the convexity of the entropy function to 
hold, the symmetric matrix w q is diagonalized by 


D = qyW g q v = q T vWv 


diag 


(R p_ p_ p_ p r \ 

\p’T'T'T' (7-I )T 2 ) ’ 


(B6) 


where 


U r — 
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V\ 
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P 
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0 

Rp ) 
7-1 / 


(B 7 ) 


A sufficient condition to ensure that C^WqQ > 0 is to ensure that all diagonal terms 
above are positive. This is satisfied by p, T > 0 . The inverse, S~^ is also calculated 
using an expansion, q w = q v v w . where 


( £- 

pu 1 
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pU3 

2 RT p— puf— pU 2 ~ pu^—2 hp 
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V 0 
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(B8) 


B.2 Viscous Stability 


To develop the viscous coefficient matrices used to define the viscous fiuxes, the 
definition of the Cartesian viscous fiuxes based on the primitive variable gradients 
is examined, 


f( y )i = (o, rii, 72 , 73, TjiVj - qi) T , z = 1,2,3. 
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This can be expressed as /^* = d i jV Xj , where 
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0 |/x 0 0 0 

C] j = 0 0 // 0 0, 

0 0 0 n 0 

\ 0 I fivi I-W 2 /uv 3 K J 
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0 n 
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0 0 
\ 0 )J,V\ 


The symmetrized coefficient matrices 
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and thus take the form, 
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(B9) 


For the viscous terms to be entropy dissipative, 


w xAjW Xj > 0 


must be satisfied. The easiest way to ensure this is to create a larger coefficient 
matrix, 


C = 


Cll 

C12 

C 13 \ 


C21 

C22 

C 23 1 , 

(BIO) 

C31 

C32 

C 33 / 
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and require 


c 


xi 


(Cx,Xx2:Xx 3 )C I Cx 2 >o, VC 


c 


X3 

To show this, the Cholesky decomposition, C = LDL T of the (15 x 15) matrix is 
calculated. The diagonal terms must non-negative, 

D= ^0, -T/i, Tfi,, Tfi, T 2 k, 0, 0, Tfi, Tfi, T 2 k, 0, 0, 0, 0, T 2 re^ . (Bll) 

In the discrete proof, the viscous terms require the additional property that each 
symmetric diagonal sub-block is also positive semi-definite. This is verified by, 


cn = LiiH n L 11 , -Du = ( 0, -Tfi,Tfi,Tfi,T k 


C 22 — L 22 D 22 L 22 , D 22 — ( 0, Tfi, ~^Tfi, Tfi, T k 


(B12) 


C33 = L 33 D 33 L 33 , D 33 — [ 0,Tfl,T fl, -Tfl,T K 


From the above, it is clear that in order for the viscous conditions to be satisfied, 


T > 0 , fi(T) > 0 , k{T) > 0 . 

Thus, recall that the contribution of the viscous terms to the semi-discrete en- 
tropy decay from 4.42 is 

W T v x; V- 1 A Xi I { [v)i w T W~}B Xi [ca]V Xc w 

i= 1 i= 1 1=1 

3 3 

- XZ v [c«] (^ w ) 

i= l 1=1 
3 3 

- (• /v ^ w ) T vv x^[cn\i (^w) . 
i= 1 t= 1 

The last term is easily shown to be always negative. Above it was shown that 
C xi&iiCxi > 0 is satisfied in each direction. Thus, all that is required is 

C T vv-%XC > 0, VC- (B13) 

Examining A13, it is clear that [ca)f is composed of convex combinations of sym- 
metric, positive semi-definite matrices, so the full matrix will also be positive semi- 
definite. The diagonal matrices 'PVT 1 are positive definite and commute with [ca\^. 
They will not effect the eigenvalues of the matrix. Thus, the requirement above 
holds and the last term will always dissipate entropy. The second-to-last term on 
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the right side of the entropy decay rate above is only slightly more difficult. It is 
first reorganized using 


3 3 / D Xl w \ 7 / V[c n ] V[c 12 ] V[ciz] \ ( \ 

^[C2l] P[c 2 2] 7^[c 2 3] , 

i=1 t=l V ^3 W / V ^[C3l] ^[C32] 7 ^[c 33 ] / \ / 

where P commutes with the viscous coefficient matrices, so it can be split into 

V = VpVP and absorbed into the gradient vectors. The gradient vectors can be 
reorganized such that the center term becomes 

3 3 

( P si w ) T (^w) = (w x )^[C](w x ) v ^, 

i= 1 t=l 

where [C] is a block diagonal matrix with blocks Ci corresponding to the viscous 
coefficients at each solution point. Each of these blocks is positive semi-definite, so 
the full matrix is positive semi-definite. Thus, 

3 3 

EE {V Xi w) T V[c ie \ (P^w) > 0, (B14) 

i= i e=i 

completing the entropy stability proof for the Navier-Stokes viscous terms. 
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